Interval Truth Model

Visualizations of Priors in the ITM

Authors

Affiliation

Matthias Kloft

University of Marburg

Björn S. Siepe

University of Marburg

Daniel W. Heck

University of Marburg

1 Latent Appraisal

1.1 True Intervals

1.1.1 Informative Beta Prior on Marginal Locations and Widths

We are mainly interested in the latent consensus \([T^{loc}_{j}, T^{wid}_{j}]^\top\). First, since we know that very wide interval widths are highly unlikely and also not meaningful in most use cases, we assign a weakly informative prior on the true interval widths on the bounded scale: \[ \begin{equation} T^{wid(0,1)}_j \sim \text{Beta}(1.2, 3) \end{equation} \] This prior has an expected value of \(.29\) and a mode of \(.09\) and therefore reflects our beliefs about the marginal true interval widths more adequately than a uniform prior.

Second, conditional on a particular interval width, we do not believe that particular interval locations are more likely than others. Therefore, we assign an uninformative prior to an auxiliary shifting weight parameter \(s_j\), which is subsequently used to compute the actual true locations on the bounded scale: \[ \begin{equation} \begin{split} s_j &\sim \text{Beta}(1, 1), \\ T^{loc(0,1)}_j &= s_j (1 - T^{wid(0,1)}_j) + \frac{T^{wid(0,1)}_j}{2}. \\ \end{split} \end{equation} \] This means that for a given interval width, we take what is left of the response scale and multiply it by the shifting weight \(s_j\), which results in the lower bound for this particular interval. To arrive at the interval location, we add half of the respective interval width to this lower bound.

Third, we transform the true interval from the bounded simplex to the unbounded bivariate scale via the isometric log-ratio function: \[ \begin{equation} \boldsymbol T^* = \text{ILR}\Bigg(\bigg[ T^{loc(0,1)}_j - \frac{T^{wid(0,1)}_j}{2}, \; T^{wid(0,1)}_j, \; T^{loc(0,1)}_j + \frac{T^{wid(0,1)}_j}{2} \bigg]^\top\Bigg). \end{equation} \]

1.1.2 Sample From Priors

# priors on marginal locations and widths
df <- data.frame(Tr_wid_splx = rbeta(1e5, 1.2, 3))
df$Tr_loc_splx <- 
  (1 - df$Tr_wid_splx) * rbeta(1e5,1,1) + df$Tr_wid_splx / 2
# transform to bivariate normal
Tr_bvn <- simplex_to_bvn(
  cbind(
    df$Tr_loc_splx - .5 * df$Tr_wid_splx,
    df$Tr_wid_splx,
    1 - (df$Tr_loc_splx + .5 * df$Tr_wid_splx)))

df <- cbind(df, Tr_loc = Tr_bvn[, 1], Tr_wid = Tr_bvn[, 2])

# transform back to simplex
Tr_splx <- bvn_to_simplex(
  cbind(df$Tr_loc, df$Tr_wid))
df$Tr_loc_splx <- Tr_splx[,1] + .5 * Tr_splx[,2]
df$Tr_wid_splx <- Tr_splx[,2]

1.1.3 Marginal Distributions

df %>%
  ggplot() +
  geom_histogram(
    aes(Tr_loc_splx),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Tr_wid_splx),
               fill = "purple",
               alpha = .3,
               binwidth = .005) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

df %>%
  ggplot() +
  geom_histogram(
    aes(Tr_loc),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Tr_wid),
               fill = "purple",
               alpha = .3,
               binwidth = .005, ) +
  scale_x_continuous(
    limits = c(NA, NA),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

1.1.4 Joint Distribution

1.1.5 Ternary

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Tr_loc_splx, y = Tr_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

1.1.5.1 Unbounded Scale

df %>%
  ggplot() +
  geom_point(
    aes(x = Tr_loc, y = Tr_wid),
    size = 1,
    alpha = .1,
    shape = 16
  ) +
  theme_pubr()

1.2 Lambda (Item Discernibility)

1.2.1 sigma_lambda

df$sigma_lambda <- exp(rnorm(1e5, log(.5), .5))
#median(lambda$sigma_lambda)

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_lambda),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.2.2 lambda

df$lambda_loc <- exp(rnorm(1e5, 0, df$sigma_lambda))
df$lambda_wid <- exp(rnorm(1e5, 0, df$sigma_lambda))

df %>%
  ggplot() +
  geom_histogram(aes(x = lambda_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3 E (Person Proficiency)

1.3.1 mu_E

Then narrow Student t prior ensures that the population proficiency is regularized to 1.

df$mu_E <- rnorm(1e5)

df %>%
  ggplot() +
  geom_histogram(aes(x = exp(mu_E)),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.2 sigma_E

df$sigma_E <- exp(rnorm(1e5, log(.5), .5))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_E),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.3 E

df$E_loc <- exp(rnorm(1e5, mean = df$mu_E, sd = df$sigma_E))
df$E_wid <- exp(rnorm(1e5, mean = df$mu_E, sd = df$sigma_E))

df %>%
  ggplot() +
  geom_histogram(aes(x = E_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.4 Lambda and E joint

df %>%
  ggplot() +
  geom_histogram(aes(x = 1 / lambda_loc / E_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "SD of Latent Appraisal", y = "Density") +
  theme_pubr()

1.4 Omega (Residual Correlation)

beta_param <- 2
df$omega <- rbeta(1e5, shape1 = beta_param, shape2 = beta_param)*2-1

df %>%
  ggplot() +
  geom_histogram(aes(x = omega),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(-1, 1), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.5 Latent Appraisals

1.5.1 Marginal Appraisals

A <- rbvnorm(1e5, 
             mean1 = df$Tr_loc, 
             mean2 = df$Tr_wid,
             sd1 = 1 / df$lambda_loc / df$E_loc,
             sd2 = 1 / df$lambda_wid / df$E_wid,
             cor = df$omega)

df$A_loc <- rnorm(1e5, mean = df$Tr_loc, sd = 1 / df$lambda_loc / df$E_loc)
df$A_wid <- rnorm(1e5, mean = df$Tr_wid, sd = 1 / df$lambda_wid / df$E_wid)

A_splx <- bvn_to_simplex(df[,c("A_loc", "A_wid")])
df$A_loc_splx <- A_splx[,1] + .5 * A_splx[,2]
df$A_wid_splx <- A_splx[,2]

df %>%
  ggplot() +
  geom_density(aes(A_loc_splx),
               fill = "lightblue",
               alpha = .3,
               shape = 16) +
  geom_density(aes(A_wid_splx),
               fill = "purple",
               alpha = .3,
               shape = 16) +
  scale_x_continuous(limits = c(0, 1), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Appraisal Locations / Widths", y = "Density") +
  theme_pubr()

1.6 Joint Appraisals

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = A_loc_splx, y = A_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

2 Biases

2.1 Shifting Bias

2.1.1 sigma_b

df$sigma_b <- exp(rnorm(1e5, log(.5), 1))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_b),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

2.1.2 b

df$b_loc <- rnorm(1e5, mean = 0, sd = df$sigma_b)
df$b_wid <- rnorm(1e5, mean = 0, sd = df$sigma_b)

df %>%
  ggplot() +
  geom_histogram(aes(x = b_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(-5, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

b_splx <- bvn_to_simplex(df[,c("b_loc", "b_wid")])
df$b_loc_splx <- b_splx[,1] + .5 * b_splx[,2]
df$b_wid_splx <- b_splx[,2]

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = b_loc_splx, y = b_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

2.2 a (Scaling Bias)

2.2.1 sigma_a

df$sigma_a <- exp(rnorm(1e5, log(.5), .5))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_a),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

2.2.2 a

df$a_loc <- exp(rnorm(1e5, mean = 0, sd = df$sigma_a))

df %>%
  ggplot() +
  geom_histogram(aes(x = a_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

3 Responses

3.1 Marginal Responses

df$Y_loc <- df$A_loc * df$a_loc + df$b_loc
df$Y_wid <- df$A_wid + df$b_wid

Y_splx <- bvn_to_simplex(df[,c("Y_loc", "Y_wid")])
df$Y_loc_splx <- Y_splx[,1] + .5 * Y_splx[,2]
df$Y_wid_splx <- Y_splx[,2]

df %>%
  ggplot() +
  geom_histogram(
    aes(Y_loc_splx),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Y_wid_splx),
               fill = "purple",
               alpha = .3,
               binwidth = .005) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

3.2 Joint Responses

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Y_loc_splx, y = Y_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

---

4 Recovery Check

4.1 Define Parameter Bounds for Variances

# compute a benchmark for the mean and SD of the parameters
mean_benchmark <- simplex_to_bvn(c(.4, .2, .4), type = "ilr")
sd_benchmark_loc <- simplex_to_bvn(c(.98, .01, .01), type = "ilr")
sd_benchmark_wid <- simplex_to_bvn(c(.495, .01, .495), type = "ilr")
# mean for Tr_loc
mu_Tr_loc <- mean_benchmark[1]
# mean for Tr_wid
mu_Tr_wid <- mean_benchmark[2]
# SD forTr_loc
sigma_Tr_loc <- sd_benchmark_loc[1] / 4
# SD Tr_wid
sigma_Tr_wid <- abs(sd_benchmark_wid[2] - mean_benchmark[2]) / 4

# SDs for other parameters
sigma_lambda_E_loc <- .3
sigma_lambda_E_wid <- .3
sigma_a_loc <- .3
sigma_b_loc <- sigma_Tr_loc / 3
sigma_b_wid <- sigma_Tr_wid / 3

4.2 Generate Data for One Replication

n_respondents <- 50
n_items <- 20

sim_data <-
  generate_itm_data_sim_study(
    n_respondents = n_respondents,
    n_items = n_items,
    mu_Tr_loc = mu_Tr_loc,
    mu_Tr_wid = mu_Tr_wid,
    sigma_Tr_loc = sigma_Tr_loc,
    sigma_Tr_wid = sigma_Tr_wid,
    sigma_lambda_E_loc = sigma_lambda_E_loc,
    sigma_lambda_E_wid = sigma_lambda_E_wid,
    sigma_a_loc = sigma_a_loc,
    sigma_b_loc = sigma_b_loc,
    sigma_b_wid = sigma_b_wid,
    omega_beta = 3
  )
responses <- sim_data$responses

4.2.1 Plot Data

samples <- gather_values(
  lower = responses$x_L,
  upper = responses$x_U,
  item_id = responses$jj,
  step_size = 0.01
)
samples %>% 
  ggplot(aes(x = samples)) +
  geom_density(fill = "purple", alpha = .5) +
  facet_wrap(~item_id) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

4.2.2 Plot True Intervals as Error Bars

data.frame(
  idx = 1:n_items,
  lower = sim_data$parameters$Tr_L,
  upper = sim_data$parameters$Tr_U
) %>%
  ggplot(aes(y = idx)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.3) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = seq(0, 1, .25),
    expand = expansion()
  ) +
  labs(x = "Latent Truth",
       y = "Item") +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

4.2.3 Plot True Intervals as Ternary Plot

data.frame(
  idx = 1:n_items,
  Tr_L = sim_data$parameters$Tr_L,
  Tr_U = sim_data$parameters$Tr_U) %>% 
  mutate(
  Tr_loc = Tr_L + (Tr_U - Tr_L) / 2,
  Tr_wid = Tr_U - Tr_L
) %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Tr_loc, y = Tr_wid),
    size = 4,
    alpha = .5,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

4.3 Fit Model

4.3.1 Stan Data

### Stan data declaration
I <- length(unique(responses$ii))
J <- length(unique(responses$jj))
N <- nrow(responses)
ii <- responses$ii
jj <- responses$jj
nn <- c(1:N)
Y_splx <- cbind(responses$x_splx_1, responses$x_splx_2, responses$x_splx_3) %>%
  as.matrix()

#Y_splx[is.na(Y_splx)]

# sanity check: all rows sum to 1
# check <- apply(Y_splx, 1, sum) %>% as.data.frame()
# table(check)

## stan data list
stan_data <- list(
  I = I,
  J = J,
  N = N,
  ii = ii,
  jj = jj,
  nn = nn,
  Y_splx = Y_splx,
  link = 1
)

4.3.2 Compile Model

# Choose model to fit
model_name <- "itm_simulation_v2_beta"
# Compile model
model <-
  cmdstanr::cmdstan_model(
    stan_file = here("src", "models", paste0(model_name, ".stan")),
    pedantic = TRUE,
    quiet = FALSE
  )

4.3.3 Run Sampler

# number of MCMC chains
n_chains <- 4
# Run sampler
fit_prior_check <- model$sample(
  data = stan_data,
  seed = 2023,
  chains = n_chains,
  parallel_chains = n_chains,
  iter_warmup = 500,
  iter_sampling = 500,
  refresh = 500,
  adapt_delta = .99,
  init = .1
)

Running MCMC with 4 parallel chains...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 228.9 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 234.9 seconds.
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 240.2 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 280.9 seconds.

All 4 chains finished successfully.
Mean chain execution time: 246.3 seconds.
Total execution time: 281.1 seconds.

# save fit_prior_check
fit_prior_check$save_object(
  file =  here("fits", paste0(model_name, "_fit_prior_check.rds")))

# load  fit
fit_prior_check <- readRDS(
  file =  here("fits", paste0(model_name, "_fit_prior_check.rds")))

---

4.3.4 Get Estimates

estimates_summary <- fit_prior_check$summary()

4.3.5 Sampler Diagnostics

fit_prior_check$diagnostic_summary()

$num_divergent
[1] 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 0.6690134 0.7145780 0.6690593 0.6188122

# sampler diagnostics
fit_prior_check$sampler_diagnostics(format = "df") %>% 
  psych::describe(quant = c(.05,.95),) %>%
  round(2) %>%  
  as.data.frame() %>% 
  dplyr::select(median, min, Q0.05, Q0.95,  max) %>% 
  .[-c(7:9),]

              median    min  Q0.05  Q0.95    max
treedepth__     7.00   6.00   6.00   7.00   7.00
divergent__     0.00   0.00   0.00   0.00   0.00
energy__      540.30 469.52 502.22 579.02 624.73
accept_stat__   0.99   0.91   0.97   1.00   1.00
stepsize__      0.04   0.04   0.04   0.05   0.05
n_leapfrog__  127.00  63.00 127.00 127.00 255.00

# convergence diagnostics
convergence_summary <- 
  fit_prior_check$draws(format = "df") %>%
  summarise_draws(.x = ., "rhat", "ess_bulk", "ess_tail") %>%
  remove_missing() %>%
  select(-variable) %>%
  psych::describe(., quant = c(.05, .95)) %>%
  as.data.frame() %>%
  select(mean, median, Q0.05, Q0.95, min, max)

convergence_summary %>% round(3)

             mean   median    Q0.05    Q0.95     min      max
rhat        1.002    1.002    1.000    1.005   0.999    1.010
ess_bulk 1956.513 1861.096  938.414 3327.645 417.098 4357.585
ess_tail 1467.217 1476.239 1161.061 1756.521 682.151 2000.306

4.3.6 Effective sample size (ESS) & Rhat Plots

# color scheme
color_scheme_set(scheme = "purple")

# Effective sample sizes
plot_neff <-
  mcmc_neff_hist(bayesplot::neff_ratio(fit_prior_check), binwidth = .01) +
  labs(title = "A") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  )
# Rhat
plot_rhat <-
  bayesplot::mcmc_rhat_hist(bayesplot::rhat(fit_prior_check)) +
  labs(title = "B") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  ) +
  yaxis_text(on = TRUE)

# Combined plot
plot_diagnostics <- gridExtra::grid.arrange(plot_neff, plot_rhat, ncol = 2)

4.4 Parameter Plots

# color scheme
color_scheme_set(scheme = "purple")

true_parameters <- sim_data$parameters

4.4.1 Person Competence: Location

plot_E_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("E_loc"), true = true_parameters$E_loc) +
  labs(
    subtitle = "Person Competence: Location",
  ) 
plot_E_loc

4.4.2 Person Competence: Width

plot_E_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("E_wid"), true = true_parameters$E_wid) +
  labs(subtitle = "Person Competence: Width")
plot_E_wid

4.4.3 Person Scaling Bias: Location

plot_a_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("a_loc"), true = true_parameters$a_loc) +
  labs(subtitle = "Person Scaling Bias: Location") +
  theme(
    text = element_text(
      family = "serif",
      size = 16,
      colour = 1
    ),
    plot.title = element_text(size = 20,
                              face = "bold"),
    plot.margin = unit(c(.1, .1, .1, .1), "cm"),
    axis.title = element_text(size = 16, color = 1),
    axis.title.x = element_text(margin = margin(t = 5)),
    axis.title.y = element_text(margin = margin(r = 5)),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x = element_text(margin = margin(t = 3))
  )
plot_a_loc

4.4.4 Person Shifting Bias: Location

plot_b_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("b_loc"), true = true_parameters$b_loc) +
  labs(subtitle = "Person Shifting Bias: Location") 
plot_b_loc

4.4.5 Person Shifting Bias: Width

plot_b_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("b_wid"), true = true_parameters$b_wid) +
  labs(subtitle = "Person Shifting Bias: Width")
plot_b_wid

4.4.6 Item Discernability: Location

plot_lambda_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("lambda_loc"),
                       true = true_parameters$lambda_loc)
plot_lambda_loc

4.4.7 Item Discernability: Width

# plot
plot_lambda_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("lambda_wid"),
                       true = true_parameters$lambda_wid) +
  labs(subtitle = "Item Discernability: Width")
plot_lambda_wid

4.4.8 Residual Correlation

plot_omega <-
  mcmc_recover_scatter(x = fit_prior_check$draws("omega"),
                       true = true_parameters$omega) +
  xlim(-1, 1) +
  ylim(-1, 1) +
  labs(subtitle = "Residual Correlation")
plot_omega

4.4.9 Latent Truth Intervals

# plot
plot_Tr_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("Tr_loc"),
                       true = true_parameters$Tr_loc) +
  labs(subtitle = "Item True Location")
plot_Tr_loc

# plot
plot_Tr_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("Tr_wid"),
                       true = true_parameters$Tr_wid) +
  labs(subtitle = "Item True Width")

plot_Tr_wid

# ITM estimates
latent_truth_est_ilr <- data.frame(
  idx = 1:J,
  type = "estimated",
  Tr_L_est = fit_prior_check$summary("Tr_splx") %>% 
    as.data.frame() %>% 
    pull(median) %>% 
    .[1:J],
  Tr_wid_est = fit_prior_check$summary("Tr_splx") %>% 
    as.data.frame() %>% 
    pull(median) %>% 
    .[(J + 1):(J * 2)]
) %>%
  mutate(Tr_loc_est = Tr_L_est + Tr_wid_est / 2, 
         Tr_U_est = Tr_L_est + Tr_wid_est)

# compute means in the unbounded space
means_ilr <-
  cbind(jj, simplex_to_bvn(Y_splx)) %>%
  as.data.frame() %>%
  dplyr::group_by(jj) %>%
  dplyr::summarise(across(everything(),mean, na.rm = TRUE)) %>%
  dplyr::ungroup() %>%
  dplyr::select(-jj) %>%
  bvn_to_simplex() %>%
  mutate(
    idx = 1:J,
    type = "simple_mean",
    Tr_L_sm = x_1,
    Tr_U_sm = 1 - x_3
  ) %>%
  select(-c(x_1, x_2, x_3))


# true parameters
latent_truth_true <- data.frame(
  idx = 1:J,
  type = "true",
  Tr_L_true = true_parameters$Tr_L,
  Tr_U_true = true_parameters$Tr_U
)

df_interval_plot_ilr <- full_join(latent_truth_est_ilr, latent_truth_true) %>%
  full_join(means_ilr) %>% 
  mutate(type = factor(type, levels = c("estimated", "true")))

4.4.9.1 Interval Plot

cols <- c("True" = "grey70","ITM (ILR)" = "red", "Mean" = "blue")

df_interval_plot_ilr %>%
  ggplot(aes(y = idx)) +
    geom_errorbarh(
    aes(xmin = Tr_L_true, xmax = Tr_U_true, col = "True"), 
    height = 0,
    linewidth = 5
  ) +
  geom_errorbarh(
    aes(xmin = Tr_L_sm, xmax = Tr_U_sm, col = "Mean"),
    height = 0,
    linewidth = 3,
    alpha = .4
  ) +
  geom_errorbarh(
    aes(xmin = Tr_L_est, xmax = Tr_U_est, col = "ITM (ILR)"),
    height = 0,
    linewidth = 2,
    alpha = .5
    ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = seq(0, 1, .25),
    expand = expansion()
  ) +
    scale_y_discrete(
    expand = expansion(add = 1)
  ) +
  scale_color_manual(values = cols) +
  labs(x = "Latent Truth", y = "Item") +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

5 Prior vs. Posterior

5.1 True Intervals

latent_truth <- data.frame(
  Tr_loc_splx_post = fit_prior_check$draws("Tr_loc_splx") %>%
    unlist() %>%
    as.vector(),
  Tr_wid_splx_post = fit_prior_check$draws("Tr_wid_splx") %>%
    unlist() %>%
    as.vector()
) %>%
  mutate(
    jj = factor(rep(1:J, each = nrow(.) / J)),
    Tr_wid_splx_prior = rbeta(nrow(.), 1.2, 3),
    Tr_loc_splx_prior = (1 - Tr_wid_splx_prior) * rbeta(nrow(.), 1, 1) + Tr_wid_splx_prior / 2
  )

Tr_bvn <- simplex_to_bvn(
  cbind(
    latent_truth$Tr_loc_splx_prior - .5 * latent_truth$Tr_wid_splx_prior,
    latent_truth$Tr_wid_splx_prior,
    1 - (
      latent_truth$Tr_loc_splx_prior + .5 * latent_truth$Tr_wid_splx_prior
    )
  )
)
latent_truth <- cbind(latent_truth,
                      Tr_loc_prior = Tr_bvn[, 1], Tr_wid_prior = Tr_bvn[, 2])


latent_truth %>%
  ggplot() +
  geom_abline(intercept = 0,
              slope = 2,
              alpha = .7) +
  geom_abline(intercept = 2,
              slope = -2,
              alpha = .7) +
  geom_density_2d_filled(
    aes(x = Tr_loc_splx_prior, 
        y = Tr_wid_splx_prior),
    binwidth = .05, 
    alpha = .3, ) +
  geom_point(
    aes(x = Tr_loc_splx_post, y = Tr_wid_splx_post, col = jj),
    size = .3,
    alpha = .3,
    shape = 16
  ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  labs(x = "Location", y = "Width") +
  # supress legends
  guides(col = FALSE, fill = FALSE) +
  theme_pubr()

5.2 Hyper-Priors

5.2.1 sigma_lambda: Location

colors <- c("prior" = "lightblue", "posterior" = "purple")

sigma_lambda_loc <- data.frame(posterior = fit_prior_check$draws("sigma_J[1]", format = "list") %>% unlist() %>%  as.vector())
sigma_lambda_loc$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_lambda_loc %>%
  ggplot() +
  geom_density(aes(prior), fill = "lightblue", alpha = 1) +
  geom_density(aes(exp(posterior * .5 + log(.5))), fill = "purple", alpha = .3) +
  scale_x_continuous(limits = c(NA, 3), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.2 sigma_lambda: Width

sigma_lambda_wid <- data.frame(
  posterior = fit_prior_check$draws("sigma_J[2]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_lambda_wid$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_lambda_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.3 sigma_E: Location

sigma_E_loc <- data.frame(
  posterior = fit_prior_check$draws("sigma_I[1]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_E_loc$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_E_loc %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.4 sigma_E: Width

sigma_E_wid <- data.frame(
  posterior = fit_prior_check$draws("sigma_I[2]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_E_wid$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_E_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.5 mu_E: Location

mu_E_loc <- data.frame(
  posterior = fit_prior_check$draws("mu_I[1]", format = "list") %>% unlist() %>%  as.vector()
)
mu_E_loc$prior <- rnorm(nrow(sigma_lambda_loc))

mu_E_loc %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(posterior),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(-4, 4),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.6 mu_E: Width

mu_E_wid <- data.frame(
  posterior = fit_prior_check$draws("mu_I[2]", format = "list") %>% unlist() %>%  as.vector()
)
mu_E_wid$prior <-  rnorm(nrow(sigma_lambda_loc))

mu_E_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(posterior),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(-4, 4),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()